El siguiente documento corresponde a un trabajo de investigación del curso Visualización de datos y análisis estadístico del Posgrado Analista de Datos de la Universidad Cenfotec.
“Estadisticas y Modelado de la distribución biogeográfica de Pumas y Jaguares de un sector del continente americano”
https://www.inaturalist.org/
https://biodiversityinformatics.amnh.org/
library(dplyr)
library(tibble)
library(tidyr)
library(sp)
library(DT)
library(leaflet)
library(leaflet.extras)
library(raster)
library(rasterVis)
library(rgdal)
library(ggplot2)
library(plotly)
library(LaCroixColoR)
library(sf)
library(here)
library(devtools)
library(proto)
library(highcharter)
library(rgeos)
library(velox)
library(usdm)
library(gtools)
library(corrplot)
library(Hmisc)
library(dismo)
En esta sección se realizan una serie de procesos que conllevan al desarrollo de las bases de datos finales
Carga de la base de datos y de las variables que la conforman
df <- read.csv("E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/wildcats_df.csv", stringsAsFactors=F, sep=",", na.strings=c("","NA"), encoding="UTF-8")
colnames(df)
## [1] "X.U.FEFF.common_name" "scientific_name" "latitude"
## [4] "longitude" "created_at" "updated_at"
## [7] "place_country_name"
Cambio de nombre de las columnas
df <- df[,c(1:5,7)]
names(df)[1] <- "Nombre_Comun"
names(df)[2] <- "Nombre_Cientifico"
names(df)[3] <- "Latitud"
names(df)[4] <- "Longitud"
names(df)[5] <- "Fecha_Registro"
names(df)[6] <- "Pais"
colnames(df)
## [1] "Nombre_Comun" "Nombre_Cientifico" "Latitud"
## [4] "Longitud" "Fecha_Registro" "Pais"
Se imprime un resumen de datos por especie de modo que se pueda obtener un conteo inicial
df%>%
dplyr::group_by(Nombre_Comun)%>%
dplyr::summarise(length(Nombre_Comun))
Se realiza una correccion de datos en la columna de “Nombre_Comun” de forma que se estandaricen algunos registros que figuran como duplicados
ndf <- df%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Norte", "Puma",Nombre_Comun))%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Sur", "Puma",Nombre_Comun))
Y por ultimo se grafica un resumen final de datos por especie
count_sp <- ndf%>%
group_by(Nombre_Comun)%>%
summarise(Count=length(Nombre_Comun))
plot_ly(data=count_sp, x = ~Nombre_Comun, y = ~Count, type = 'bar', text = ~Count, textposition = 'auto',
marker = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
line = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
width = 1.5))) %>%
layout(title = "Cantidad de especies registrados en el dataframe inicial",
xaxis = list(title = "Especies"),
yaxis = list(title = "Cantidad"))
En el siguiente mapa se imprime, en la figura de puntos, la totalidad de los datos
cf <- colorFactor(c("#ffa500", "#13ED3F"), domain=ndf$Nombre_Comun)
m <- leaflet(df) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=ndf$Nombre_Comun, group = ndf$Nombre_Comun ,weight = 3, radius=4,
color=cf(ndf$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = ndf$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m
El area de estudio con la que se pretende trabajar corresponde a Centro y Sudamrica, asi como el Caribe.
Algunas de las coberturas que se utilizan para conocer el habitat de las especies son variables climáticas que se derivan de los datos proporcionados por el Panel Intergubernamental sobre Cambio Climático y se produjeron utilizando Interpolación de basada en lecturas tomadas en estaciones meteorológicas de todo el mundo desde 1961 hasta 1990.
Los coberturas corresponden a:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*ecoreg : Ecoregiones
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_ann : precipitacion, anual
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*tmn6190_ann : temperatura promedio, anual
*tmp6190_ann : temperatura minima, anual
*tmx6190_ann : temperatura maxima, anual
*vap6190_ann : presion de vapor, anual
r_files <- list.files(here::here("Datasets", "coverages"),
full.names = T)
r_files
## [1] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/cld6190_ann.asc"
## [2] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/dtr6190_ann.asc"
## [3] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/ecoreg.asc"
## [4] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/frs6190_ann.asc"
## [5] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/h_dem.asc"
## [6] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_ann.asc"
## [7] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l1.asc"
## [8] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l10.asc"
## [9] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l4.asc"
## [10] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l7.asc"
## [11] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmn6190_ann.asc"
## [12] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmp6190_ann.asc"
## [13] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmx6190_ann.asc"
## [14] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/vap6190_ann.asc"
Se grafican los raster que representan las variables ambientales con las que se utilizaran en el modelado de distribucion
st <- raster::stack(r_files)
raster::plot(st)
La variable ejemplicada, ecoregion, posee las siguientes categorias:
1: Mediterranean Forests, Woodlands & Scrub 2: Tropical & Subtropical Dry Broadleaf Forests 3: Temperate Broadleaf & Mixed Forests 4: Tropical & Subtropical Coniferous Forests 5: Temperate Grasslands, Savannas & Shrublands 6: Mangroves 7: ( ) 8: Montane Grasslands & Shrublands 9: Tropical & Subtropical Grasslands, Savannas & Shrublands 10: Tropical & Subtropical Moist Broadleaf Forests 11: Flooded Grasslands & Savannas 12: Montane Grasslands & Shrublands 13: Magellanic subpolar forests 14: Rock and Ice 15: ( )
ecoreg = raster(r_files[3])
raster::plot(ecoreg, main="Ecoregiones del area de estudio")
raster::hist(ecoreg, main="Histograma de las Ecoregiones del area de estudio", col="green")
A continuacion, se carga el un archivo espacial de tipo vectorial correspondiente a los paises que conforman el area de estudio.
area_estudio <- sf::st_read(dsn="E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/map.geojson", layer="map", quiet=TRUE)
En el siguiente mapa se muestra los paises que conforman el area de estudio
map.area <- plot_geo(area_estudio, name=~pais, color = I("#89C5DA"))%>% hide_legend()%>%
layout(title = "Mapa paises que conforman el area de estudio")
map.area
De previo se realiza una selección de la base de datos de las especies respecto del área de estudio
coords <- data.frame(x=ndf$Longitud, y=ndf$Latitud)
prj <- CRS("+proj=longlat +datum=WGS84 +no_defs")
points <- SpatialPointsDataFrame(coords, ndf, proj4string = prj)
wildcats <- crop(points,ecoreg)
wildcats <- data.frame(wildcats)
wildcats <- wildcats[,1:6]
Para luego visualizar en un mapa dinamico los puntos de observacion de las especies a estudiar
cf2 <- colorFactor(c("#ffa500", "#13ED3F"), domain=wildcats$Nombre_Comun)
m2 <- leaflet(wildcats) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=wildcats$Nombre_Comun, group = wildcats$Nombre_Comun ,weight = 3, radius=4,
color=cf2(wildcats$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = wildcats$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m2
Primero se realiza una inclusion de nuevas variables relacionadas con el tema de temporalidad
wildcats$Hora_Registro <- strftime(wildcats$Fecha_Registro, format="%H") #col7
wildcats$Fecha_Registro <- as.POSIXct(wildcats$Fecha_Registro, tz="", format="%Y-%m-%d") #col5
wildcats$Y.M.R <- strftime(wildcats$Fecha_Registro, format="%Y-%m") #col8
wildcats$M.R <- format(wildcats$Fecha_Registro, "%m") #col9
wildcats$Y.R <- format(wildcats$Fecha_Registro, "%Y") #col10
wildcats <- wildcats[, colnames(wildcats)[c(1:4,6,5,8,10,9,7)]] #Reorganizacion de las columnas
y posteriormente se procede a visualizar la base de datos final
datatable(wildcats, extensions = c('Buttons','FixedColumns', 'RowReorder'),
filter = list(position = 'top', clear = FALSE),
options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), dom = 't',
scrollX = TRUE,
fixedColumns = TRUE,
rowReorder = TRUE, order = list(c(0 , 'asc'))))
En este primer apartado se realizan analisis sobre estadísticas espacio-temporales respecto de los pumas y jaguares de un sector del continente americano
En los siguientes graficos se visualizan cantidad de especies por pais
En esta primera grafica de este subapartado, se muestra la cantidad por pais en relacion con el mes en que las observaciones ocurrieron
p_m <- plot_ly(data=wildcats, x= ~Pais, split= ~Nombre_Comun, color= I(cf2(wildcats$Nombre_Comun)),
frame=wildcats$M.R, type = "histogram", alpha = 0.6)%>%
layout(barmode = "overlay", title = "Cant de especies por pais y mes de observacion")
p_m%>%animation_opts(frame= 2000, redraw = TRUE)%>%
animation_slider(currentvalue = list(prefix = "Mes ", font = list(color="red")))
Para poder representar la cantidad de especies por pais se realiza un preseleccion del dato por especie
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())
Y posteriormente se procede a graficar la informacion primero de jaguares
data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Jaguares") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.jaguars, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas=FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE)
Y posteriormente de pumas
data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Pumas") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.pumas, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas= FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE)
En este subapartado, se analizan cantidades de las especies observados en relacion a la temporalida de las observaciones
Al igual como se hizo en un punto anterior, para poder representar la cantidad de especies segun periodo de observacion, se realiza un preseleccion del dato por especie
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
ct.jaguars <- ts(data=c.jaguars, start = c(2011,10), end = c(2019,12), frequency = 12)
ct.pumas <- ts(data=c.pumas, start = c(2011,10), end = c(2019,12), frequency = 12)
Y se grafica posteriormente una serie de tiempo en una linea de cantidades de especies observadas
highchart(type = "stock") %>%
hc_title(text = "Cantidad de especies") %>%
hc_subtitle(text = "en la linea de tiempo de observacion") %>%
hc_add_series(ct.jaguars[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de jaguares observados a la fecha", color = "#ffa500")%>%
hc_add_series(ct.pumas[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de pumas observados a a fecha", color = "#13ED3F")%>%
hc_xAxis(title= list(text='Fechas'), type='datetime')%>%
hc_yAxis(title= list(text='Cantidad de especies observadas'), opposite = FALSE)
Asi mismo, se realiza una descomposicion de la serie de tiempo para observar el comportamiento de la linea temporal segun la observada, la tendencia, la estacional y la aleatoria.
decomp_jaguars <- decompose(ct.jaguars[,2])
decomp_pumas <- decompose(ct.pumas[,2])
plot(decomp_jaguars)
par(adj = 1)
title(sub = "Jaguares", cex.sub = 1.5, font.sub = 3, col.sub = "#ffa500")
plot(decomp_pumas, sub="Pumas")
par(adj = 1)
title(sub = "Pumas", cex.sub = 1.5, font.sub = 3, col.sub = "#13ED3F")
En este apartado se analizan la ubicacion de las especies observadas en terminos distribucion, densidad por espacio ocupado (agrupamiento), ademas de su ubicacion con relacion a la ecoregion que ocupan.
En estos primeros mapas se visualiza la correspondiente ubicacion de los individuos de las diferentes especies, asi como la dinamica de emplazamiento en relacion al mes en que se observaron
e.m <- ggplot(wildcats, aes(x=Longitud, y=Latitud, frame=M.R, text = paste("Pais:", Pais))) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 1)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Distribucion de las especies observadas segun mes")
ggplotly(e.m)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
Al igual que los anteriores mapas de este subapartado, se se visualiza la correspondiente ubicacion de los individuos de las diferentes especies, haciendo referencia al patron de distribucion eliptico que estos individuos siguen en la dinamica de emplazamiento respecto del mes en que se observaron
pd <- qplot(data = wildcats, x=Longitud, y=Latitud) +
stat_ellipse(geom = "polygon", alpha = 1/2, aes(fill = Nombre_Comun, frame=M.R))+
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Patron de Distribucion de las Especies")
ggplotly(pd)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
En este par de mapas sobre cada especie estudiada, se visualiza el nivel de densidad que muestra el emplazamiento de los individuos en los diferentes meses en que estos fueron observaron
dn <- ggplot(wildcats, aes(x=Longitud, y=Latitud)) +
borders("world", regions=".", xlim = c(-120, -30), ylim = c(-50, 25))+
stat_density2d(aes(fill = stat(level), frame=wildcats$M.R), geom="polygon", alpha=0.3) +
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2, alpha=0.3)+
facet_grid(cols = vars(Nombre_Comun)) +
scale_fill_viridis_c(option = "viridis") +
theme(legend.position = "magma") +
labs(title = "Densidad segun distribucion por especie")
ggplotly(dn)%>%
animation_opts(
4000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
La categorias de Ecoregiones/Biomas que se observan son:
1: Mediterranean Forests, Woodlands & Scrub 2: Tropical & Subtropical Dry Broadleaf Forests 3: Temperate Broadleaf & Mixed Forests 4: Tropical & Subtropical Coniferous Forests 5: Temperate Grasslands, Savannas & Shrublands 6: Mangroves 7: ( ) 8: Montane Grasslands & Shrublands 9: Tropical & Subtropical Grasslands, Savannas & Shrublands 10: Tropical & Subtropical Moist Broadleaf Forests 11: Flooded Grasslands & Savannas 12: Montane Grasslands & Shrublands 13: Magellanic subpolar forests 14: Rock and Ice 15: ( )
En este punto se muestran tres graficas diferentes en las que se representa: 1) la ubicacion espacial de los individuos segun la ecoregion (en color naranja: jaguares, verde: pumas), 2) la ubicacion espacial de los individuos por especie diferenciados por color segun ecoregion, y 3) la cantidad de individuos por especie segun ecoregion.
ecoreg_spec <- raster::extract(ecoreg, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
ecoreg_spec <- na.omit(ecoreg_spec)
ecoreg_spec$Ecoreg <- as.character(ecoreg_spec$.)
r_eco <- plot(ecoreg)+
points(x=wildcats$Longitud, y=wildcats$Latitud, pch=20, col= cf2(wildcats$Nombre_Comun))
ggplot(ecoreg_spec, aes(x=Longitud, y=Latitud, color = ecoreg_spec$Ecoreg)) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(ecoreg_spec), shape = 21, size = 1)+
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Especies observadas segun ecoregion en la que se encuentran")
ggplot(ecoreg_spec, aes(y=Ecoreg, x=Nombre_Comun, colour = Ecoreg)) +
geom_count(alpha=0.5) +
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
labs(title = "Cantidad de especies observadas por ecoregion",
x = "Especie",
y = "Ecoregiones",
size = "")
Este analisis se realiza con el fin de evaluar cuales son las variables mas adecuadas para una modelacion espacial.
Para lograr este analisis se procede a unir de los puntos de observacion de las especies junto con los raster de las variables ambientales. En este proceso se excluye ecoregion, la cual es una variable categorica. De la base de datos wildcats unicamente se trabaja con algunas cuantas variables que tiene que ver con la especie y su ubicacion. Una vez se obtiene el resultado de la union, se procede a eliminar los valores nulos, los cuales serian aquellos pixel de las variables raster que no logran punto de union espacial con los individuos de las especies.
#Se extraen temas relevantes de las variables ambientales
select_r = dropLayer(st,3)
jp_swd <- raster::extract(select_r, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
#Se elimina las filas con valores nulos
jp_swd <- na.omit(jp_swd)
Con el anterior resultado, se procede a aplicar un analisis de correlacion multiple y de igual forma graficar el resultado de dicho analisis. En la grafica se observan una representacion con circulos que poseen distintos tamaños y colores. La escala de colores esta en funcion de valores que oscilan entre 1 y -1, donde -1 es una correlacion inversa multiproporcional, es decir, que mientras una variable aumenta la otra disminuye, a su vez, valores cercanos a 1 muestran una correlacion directa multiproporcional, es decir, que hay una fuerte correlacion en cuyo caso mientras una variable aumenta la otra variable tambien lo hace, y por su parte, valores cercanos a 0 poseen una correlacion bastante defil, es decir que no hay asociacion entre las dos variables. Con el tamaño del circulo se visualiza la magnitud de la correlacion, osea, entre más grande el círculo más cercano a -1 o 1 según sea.
mt <- cor(jp_swd[,11:23])
corrplot(mt, method = 'circle', type = 'upper')
Para realizar este el analisis se emplea la regresion lineal VIF, la cual es el factor de inflación de la varianza. Este cuantifica la intensidad de la multicolinealidad en un análisis de regresión normal de mínimos cuadrados. Proporciona un índice que mide hasta qué punto la varianza (el cuadrado de la desviación estándar estimada) de un coeficiente de regresión estimado se incrementa a causa de la colinealidad. Estamos considerando colinealidad en el contexto de un modelo estadístico que se utiliza para estimar la relación entre una variable de respuesta y un conjunto de variables predictoras (“independientes” o “explicativas”), osea, describe la situación en la que dos o más variables de predicción en un modelo estadístico están relacionadas linealmente. Dado el caso de que los factores VIF mayores que 10 implican problemas graves de multicolinealidad, el argumento th, el cual hace referencia al umbral de correlación, se establece como 10, de este modo se sabra cuales variables tienen problemas de correlacion y estas se descartaran para el modelado de distribucion.
# Analisis VIF
vif.res <- vif(x = jp_swd[,11:ncol(jp_swd)])
vif.step <- vifstep(x = jp_swd[,11:ncol(jp_swd)], th = 10)
vrs <- vif.step@results$Variables %>% as.character()
vif.step
## 4 variables from the 13 input variables have collinearity problem:
##
## pre6190_ann tmp6190_ann tmn6190_ann tmx6190_ann
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( pre6190_l10 ~ pre6190_l1 ): -0.009063609
## max correlation ( vap6190_ann ~ frs6190_ann ): -0.8508242
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 cld6190_ann 3.025670
## 2 dtr6190_ann 1.678499
## 3 frs6190_ann 4.102199
## 4 h_dem 2.447628
## 5 pre6190_l1 1.755003
## 6 pre6190_l10 5.388244
## 7 pre6190_l4 4.074227
## 8 pre6190_l7 2.896292
## 9 vap6190_ann 5.176151
vrs
## [1] "cld6190_ann" "dtr6190_ann" "frs6190_ann" "h_dem" "pre6190_l1"
## [6] "pre6190_l10" "pre6190_l4" "pre6190_l7" "vap6190_ann"
A partir del paso anterior, fue posible conocer cuales variables tienen problemas de colinealidad a la vez que se conocia cuales variables se relacionan mejor entre si, los cuales se imprimen en l tabla siguiente junto con los individuos de las especies y sus respectivas ubicaciones espaciales
jp_fnl <- jp_swd %>%
as_tibble %>%
dplyr::select(Nombre_Cientifico:Longitud, vrs)
head(jp_fnl)
Con respecto al resultado de la correlacion se grafican las relaciones entre variables con multicolinealidad:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*vap6190_ann : presion de vapor, anual)
r_corr= st[[c(1,2,4,5,7,8,9,10,14)]]
raster::pairs(r_corr)
El primer paso en este proceso de modelado es extraer los valores de los predictores en las ubicaciones de los puntos. Los predictores hacen referencia a los raster que conforman las variables ambientales que se emplearan en el modelado, los cuales son los que presentan multicolinealidad. Ademas, se crean 500 puntos aleatorios con base en los predictores, los cuales son empleados para evaluar la presencia (1) y ausencia (0) de individuos de las especies estudiadas en relacion a las variables ambientales (predictores). En este punto, se imprimen tablas para ver los resultados de presencia y ausencia, asi como un resumen general de este analisis.
pvals <- raster::extract(r_corr, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
presvals <- pvals[,11:19]
predictors <- r_corr
backgr <- randomPoints(predictors, 500)
absvals <- raster::extract(predictors, backgr)
pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals)))
head(sdmdata)
tail(sdmdata)
summary(sdmdata)
## pb cld6190_ann dtr6190_ann frs6190_ann
## Min. :0.0000 Min. :32.00 Min. : 46.0 Min. : 0.00
## 1st Qu.:0.0000 1st Qu.:58.00 1st Qu.: 99.0 1st Qu.: 0.00
## Median :1.0000 Median :63.00 Median :115.0 Median : 1.00
## Mean :0.6736 Mean :63.31 Mean :111.8 Mean : 12.22
## 3rd Qu.:1.0000 3rd Qu.:68.00 3rd Qu.:123.0 3rd Qu.: 3.00
## Max. :1.0000 Max. :83.00 Max. :166.0 Max. :198.00
## NA's :9 NA's :9 NA's :9
## h_dem pre6190_l1 pre6190_l10 pre6190_l4
## Min. : 0 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 107 1st Qu.: 22.00 1st Qu.: 28.00 1st Qu.: 22.00
## Median : 192 Median : 50.00 Median : 41.00 Median : 34.00
## Mean : 484 Mean : 51.98 Mean : 53.75 Mean : 44.43
## 3rd Qu.: 509 3rd Qu.: 75.00 3rd Qu.: 78.50 3rd Qu.: 61.00
## Max. :5455 Max. :173.00 Max. :210.00 Max. :178.00
## NA's :65 NA's :9 NA's :9 NA's :9
## pre6190_l7 vap6190_ann
## Min. : 0.00 Min. : 1
## 1st Qu.: 6.00 1st Qu.:202
## Median : 29.00 Median :255
## Mean : 42.31 Mean :226
## 3rd Qu.: 70.00 3rd Qu.:265
## Max. :194.00 Max. :308
## NA's :9 NA's :9
En este proceso los métodos toman una fórmula que identifica las variables dependientes e independientes, acompañada de un marco de datos que contienen estas variables.
Para ello se utiliza la función (glm) la cual permite ajustar modelos lineales generalizados, especificados mediante una descripción simbólica del predictor lineal y una descripción de la distribución de errores.
En un primer modelo se trabajara nuevamente con las variables ambientales con multicolinealidad:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*frs6190_ann : frecuencia de heladas, anual
*h_dem : Modelo Digital de Elevacion
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*vap6190_ann : presion de vapor, anual
m1 <- glm ( pb ~ cld6190_ann + dtr6190_ann + frs6190_ann + h_dem + pre6190_l1 + pre6190_l10 + pre6190_l4 + pre6190_l7 + vap6190_ann, data = sdmdata )
class ( m1 )
## [1] "glm" "lm"
summary ( m1 )
##
## Call:
## glm(formula = pb ~ cld6190_ann + dtr6190_ann + frs6190_ann +
## h_dem + pre6190_l1 + pre6190_l10 + pre6190_l4 + pre6190_l7 +
## vap6190_ann, data = sdmdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0442 -0.3906 0.1907 0.2840 0.9945
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.788e-01 1.527e-01 -6.409 1.98e-10 ***
## cld6190_ann 9.063e-03 1.813e-03 4.999 6.45e-07 ***
## dtr6190_ann 4.767e-03 7.515e-04 6.344 2.99e-10 ***
## frs6190_ann 3.241e-03 6.658e-04 4.867 1.25e-06 ***
## h_dem -6.392e-05 2.226e-05 -2.872 0.00414 **
## pre6190_l1 1.342e-03 4.413e-04 3.042 0.00239 **
## pre6190_l10 4.757e-03 5.166e-04 9.208 < 2e-16 ***
## pre6190_l4 -5.991e-03 5.696e-04 -10.517 < 2e-16 ***
## pre6190_l7 1.224e-03 4.280e-04 2.859 0.00431 **
## vap6190_ann 1.916e-03 3.812e-04 5.025 5.65e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1769934)
##
## Null deviance: 316.67 on 1464 degrees of freedom
## Residual deviance: 257.53 on 1455 degrees of freedom
## (67 observations deleted due to missingness)
## AIC: 1632.6
##
## Number of Fisher Scoring iterations: 2
m1
##
## Call: glm(formula = pb ~ cld6190_ann + dtr6190_ann + frs6190_ann +
## h_dem + pre6190_l1 + pre6190_l10 + pre6190_l4 + pre6190_l7 +
## vap6190_ann, data = sdmdata)
##
## Coefficients:
## (Intercept) cld6190_ann dtr6190_ann frs6190_ann h_dem pre6190_l1
## -9.788e-01 9.063e-03 4.767e-03 3.241e-03 -6.392e-05 1.342e-03
## pre6190_l10 pre6190_l4 pre6190_l7 vap6190_ann
## 4.757e-03 -5.991e-03 1.224e-03 1.916e-03
##
## Degrees of Freedom: 1464 Total (i.e. Null); 1455 Residual
## (67 observations deleted due to missingness)
## Null Deviance: 316.7
## Residual Deviance: 257.5 AIC: 1633
El algoritmo BIOCLIM calcula la similitud de una ubicación comparando los valores de las variables ambientales en cualquier ubicación con una distribución porcentual de los valores en ubicaciones conocidas de ocurrencia. Cuanto más cerca del percentil 50 (la mediana) más adecuada es la ubicación. Las colas de la distribución no se distinguen, es decir, el percentil 10 se trata como equivalente al percentil 90.
Para este caso solo se usa datos de presencia, por lo que usamos ‘presvals’ en lugar de ‘sdmdata’.
bc <- bioclim ( presvals [, c ( 'cld6190_ann', 'dtr6190_ann', 'frs6190_ann', 'h_dem', 'pre6190_l1', 'pre6190_l10', 'pre6190_l4', 'pre6190_l7', 'vap6190_ann' )])
class ( bc )
## [1] "Bioclim"
## attr(,"package")
## [1] "dismo"
bc
## class : Bioclim
##
## variables: cld6190_ann dtr6190_ann frs6190_ann h_dem pre6190_l1 pre6190_l10 pre6190_l4 pre6190_l7 vap6190_ann
##
##
## presence points: 1002
## cld6190_ann dtr6190_ann frs6190_ann h_dem pre6190_l1 pre6190_l10 pre6190_l4
## 1 66 118 0 107 72 30 34
## 2 57 99 0 31 113 110 88
## 3 60 121 1 250 14 45 14
## 5 63 123 10 757 86 43 32
## 7 56 84 1 51 48 98 25
## 9 60 119 1 275 15 46 14
## 10 58 97 0 2 121 111 87
## 12 60 95 0 2 23 58 12
## 14 66 118 0 110 72 30 34
## 15 53 116 1 200 47 130 43
## pre6190_l7 vap6190_ann
## 1 6 260
## 2 150 269
## 3 42 266
## 5 6 214
## 7 134 268
## 9 43 266
## 10 165 278
## 12 47 281
## 14 6 260
## 15 171 238
## (... ... ...)
pairs ( bc )
Todos estos objetos “modelo”, independientemente de su clase exacta, pueden utilizarse con la función de “predict” con el fin de hacer predicciones para cualquier combinación de valores de las variables independientes.
A continuación, se realizan predicciones con el objeto del modelo glm ‘m1’ y para el modelo bioclim ‘bc’, para tres registros con valores para las variables: cld6190_ann, dtr6190_ann, frs6190_ann, h_dem, pre6190_l1, pre6190_l10, pre6190_l4, pre6190_l7, vap6190_ann. Esto ultimo, se muestran en un primer dataframe.
En los resultados de la prediccion se encuentran categorias donde el 1 indica menos frecuencia para encontrar especies, 2 frecuencia media y 3 donde seria mas frecuente encontrar especies.
cld6190_ann = c(40, 60, 80)
dtr6190_ann = c(40, 100, 140)
frs6190_ann = c(50, 125, 200)
h_dem = c(1000, 3000, 6000)
pre6190_l1 = c(50, 100, 150)
pre6190_l10 = c(50, 100, 150)
pre6190_l4 = c(50, 100, 150)
pre6190_l7 = c(50, 100, 150)
vap6190_ann = c(50, 100, 250)
pdt = data.frame ( cbind (cld6190_ann, dtr6190_ann, frs6190_ann, h_dem, pre6190_l1, pre6190_l10, pre6190_l4, pre6190_l7, vap6190_ann))
pdt
predict (m1 , pdt)
## 1 2 3
## -0.1651051 0.5798087 1.3570168
predict (bc , pdt)
Hacer tales predicciones para algunos entornos puede ser muy útil para explorar y comprender las predicciones del modelo, por tanto en este punto, se utiliza la función “response” para crea gráficos de respuesta para cada variable, con las otras variables en su valor medio, en funcion del modelo bioclim.
response (bc)
Por ultimo, se representan los resultados de la prediccion en un mapa para poder visualizar la posible distribucion de especies en funcion de las variables ambientales y la presencia de individuos observados de las especies estudiadas.
En el espectro de colores del mapa se observan los valores mas altos como aquellos donde es idoneo encontrar estas especies, y al contrario en el espectro, donde es menos idoneo presencia de estas.
p <- predict (predictors,m1)
plot (p)